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^ Abstract 

R. Hirota and K. Kimura discovered integrable discretizations of the Euler and 
the La grange tops, given by birational maps. Their method is a specialization to 
the integrable context of a general discretization scheme introduced by W. Kahan 
and applicable to any vector field with a quadratic dependence on phase variables. 
According to a proposal by T. Ratiu, discretizations of the Hirota-Kimura type 
can be considered for numerous integrable systems of classical mechanics. Due to a 
^Zh remarkable and not well understood mechanism, such discretizations seem to inherit 

, ^ , the integrability for all algebraically completely integrable systems. We introduce 

an experimental method for a rigorous study of integrability of such discretizations. 
Application of this method to the Hirota-Kimura type discretization of the Clebsch 
l/-) system leads to the discovery of four functionally independent integrals of motion 

of this discrete time system, which turn out to be much more complicated than 
the integrals of the continuous time system. Further, we prove that every orbit 
of the discrete time Clebsch system lies in an intersection of four quadrics in the 
six-dimensional phase space. Analogous results hold for the Hirota-Kimura type 
discretizations for all commuting flows of the Clebsch system, as well as for the 
so(4) Euler top. 

> 

•l-H 

jh 1 Introduction 



The discretization method studied in this paper seems to be introduced in the geometric 
integration literature by W. Kahan in the unpublished notes [Kahan 1993]. It is applicable 
to any system of ordinary differential equations for x : R — > IR n with a quadratic vector 
field: 

x = Q( x ) + Bx + c, (1) 
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where each component of Q : K n —>■ K n is a quadratic form, while B e Mat nxn and 
c G W 1 . Kahan's discretization reads as 

X -^- = Q(x,Z) + -B{x + x) + c, (2) 

where 

Q(x, x) = - (q(x + x) - Q(x) - Q(x)) 

is the symmetric bilinear form corresponding to the quadratic form Q. Here and below 
we use the following notational convention which will allow us to omit a lot of indices: for 
a sequence x : Z — > K we write x for x*; and x for Xk+i- Eq. (|2j) is linear with respect to x 
and therefore defines a rational map x = f(x, e). Clearly, this map approximates the time- 
(2e)-shift along the solutions of the original differential system, so that xj~ m x(2ke). (We 
have chosen a slightly unusual notation 2e for the time step, in order to avoid appearance of 
various powers of 2 in numerous formulas; a more standard choice would lead to changing 
e i — > e/2 everywhere.) Since eq. ^ remains invariant under the interchange x <-> x with 
the simultaneous sign inversion e i— > — e, one has the reversibility property 

r 1 (x,6) = /(x,-e). (3) 

In particular, the map / is birational. 

W. Kahan applied this discretization scheme to the famous Lotka-Volterra system and 
showed that in this case it possesses a very remarkable non-spiralling property. We will 
briefly discuss this example in Sect. [2] Some further applications of this discretization 
have been explored in [Kahan and Li 1997]. 

The next, even more intriguing appearance of this discretization was in the two papers 
by R. Hirota and K. Kimura who (being apparently unaware of the work by Kahan) 
applied it to two famous integrable system of classical mechanics, the Euler top and the 
Lagrange top [Hirota and Kimura 2000, Kimura and Hirota 2000]. For the purposes of the 
present text, integrability of a dynamical system is synonymous with the existence of a 
sufficient number of functionally independent conserved quantities, or integrals of motion, 
that is, functions constant along the orbits. We leave aside other aspects of the multi-facet 
notion of integrability, such as Hamiltonian ones or explicit solution. Surprisingly, the 
Kahan-Hirota-Kimura discretization scheme produced in both the Euler and the Lagrange 
cases of the rigid body motion integrable maps. Even more surprisingly, the mechanism 
which assures integrability in these two cases seems to be rather different from the majority 
of examples known in the area of integrable discretizations, and, more generally, integrable 
maps, cf. [Suris 2003]. The case of the discrete time Euler top is relatively simple, and 
the proof of its integrability given in [Hirota and Kimura 2000] is rather straightforward 
and easy to verify by hands. As it often happens, no explanation was given in [Hirota and 
Kimura 2000] about how this result has been discovered. The "derivation" of integrals of 
motion for the discrete time Lagrange top in [Kimura and Hirota 2000] is rather cryptic 
and almost uncomprehensible. 

The present paper aims at clarifying the Hirota-Kimura integrability mechanism and 
at its application to further integrable systems. We use the term "Hirota-Kimura type 
discretization" for the Kahan's discretization in the context of integrable systems. In Sect. 
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[3] we propose a formalization of the Hirota-Kimura mechanism from [Kimura and Hirota 
2000], which will hopefully unveil its main idea and contribute towards demystifying at 
least some of its aspects. We introduce a notion of a "Hirota-Kimura basis" for a given 
map /. Such a basis $ is a set of simple (often monomial) functions, $ = (<pi, . . . , <pi), 
such that for every orbit {f l {x)} of the map / there is a certain linear combination 
C\ip\ + . . . + ci<fi of functions from $ vanishing on this orbit. As explained in Sect. |3j 
this is a new mathematical notion, not reducible to that of integrals of motion, although 
closely related to the latter. In Sect. [4] we lay a theoretical fundament for the search 
for Hirota-Kimura bases for a given discrete time system, and give a number of practical 
recipes and tricks for doing this. 

We dare to claim that the results of [Hirota and Kimura 2000] concerning the discrete 
time Euler top were originally discovered using the mechanism of Hirota-Kimura bases, 
and we present in Sect. [5] an attempt to reconstruct the way this discovery has been 
made. Sect. [6] contains the main results of this paper, namely the proof of integrability 
of the Hirota-Kimura type discretization for a further famous integrable system of the 
classical mechanics, namely for the Clebsch case of the motion of a rigid body in an ideal 
fluid. 

Our investigations are based mainly on computer experiments, which are used both for 
discovery of new results and for their rigorous proof. A search for Hirota-Kimura bases 
can be done with the help of numerical experiments based on the recipe (N) formulated 
in Sect. |4j which has a theoretical justification in Theorem [6j If the search has been 
successful and a certain set of functions $ has been identified as a Hirota-Kimura basis 
for a given map /, then numerical experiments can provide a very convincing evidence 
in favor of such a statement. A rigorous proof of such a statement turns out to be much 
more demanding. At present, we are not in possession of any theoretical proof strategies 
and are forced to verify the corresponding statements by means of symbolic computations. 
However, direct and simple-minded symbolic computations turn out to be non-feasible due 
to complexity issues. As detailed in Sect. [6| the sheer size of explicit expressions for the 
second iterate f 2 of the discrete time Clebsch system precludes symbolic manipulations, 
like solution of linear systems, as soon as those involve f 2 . Therefore our main effort 
has been put into finding the strategy of a complete and rigorous symbolical proof which 
would avoid using f 2 and would stay within the memory and performance restrictions of 
the available software and hardware. The resulting proofs are computer assisted and are 
based on symbolic computations with MAPLE [MAPLE], SINGULAR [SINGULAR] and 
FORM [FORM]. 

Our work was stimulated by a talk by T. Ratiu at the Oberwolfach Workshop "Geo- 
metric Integration" [Ratiu 2006], where an extension of the Hirota-Kimura approach to 
the Clebsch system and to the Kovalevski top has been proposed. However, no valid 
derivation of integrals was presented in T. Ratiu's talk, so that the question on the inte- 
grability of these discretizations remained open. Our work answers this question in the 
affirmative for the Clebsch system (actually, even for a whole family of Hamiltonian flows 
generated by commuting integrals of the Clebsch system). In the concluding Sect. [7j we 
discuss further perspectives of this approach and formulate a general conjecture about 
the integrability of the Hirota-Kimura type discretizations.. 
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2 Kahan's discretization of the Lotka-Volterra sys- 
tem 



As already mentioned in Sect. [TJ W. Kahan applied his general discretization scheme to 
the famous Lotka-Volterra system modelling the interaction of the predator and the prey 
populations: 

x = x(l-y), y = y(x-l). (4) 

Solutions of this system lie on closed curves in (the first quadrant of) the phase plane M 2 , 
because of the presence of the integral (conserved quantity) 

H{x,y) = x + y-\og(xy). 

Actually, system Q is Hamiltonian with respect to the Poisson bracket 

{x, y} = xy, (5) 

with the Hamilton function H: 

dH dH 
x = -xy— , y = xy— . 
ay ox 

The majority of the conventional discretization schemes produce, when applied to Q, 
spiralling solutions. Compared with solutions of the original system, this is a qualitatively 
different behavior, cf. Fig. [2] (left). The discretization proposed by Kahan reads: 



(x- x)/e = (x + x) - (xy+xy), (y - y)/e = (xy + xy) - (y + y), (6) 

Eq. (J6]) can be written as a linear system for (x, y) , 

(\ — e + ey ex \ /x\ _ f(l + e)x\ 
V -ey 1 + e-ex) \y) = {(l-e)y) ' 

which can be immediately solved, thus yielding an explicit map (x, y) = f(x,y,e): 

+ e) 2 -e(l + e)x-e(l-e)y 



x = x 



1 - e 2 - e(l - e)x + e(l + e)y ' 
(1 - e) 2 + e(l + e)x + e(l - e)y 

V 1 - e 2 - e(l - e)x + e(l + e)y ' 



(7) 



A remarkable property of the Kahan's discretization is that it apparently does not suffer 
from spiralling, solutions seem to fill out closed curves in the phase plane, cf. Fig. [2] 
(right). A (partial) explanation of this behavior was given in [Sanz-Serna 1994], where 
it was shown that the map / is Poisson with respect to the invariant Poisson bracket 
(|5| of the system Q. It is unknown whether the map [7] possesses an integral of motion, 
thus forcing all orbits to lie on smooth closed curves, as suggested by Fig. [2] (right). 
Some numerical experiments, via a deep zoom-in into certain domains of the phase plane, 
indicate that the map might be non-integrable, but a rigorous proof of a non-existence 
statement seems to be rather difficult. It might be possible with the use of technology 
described in [Gelfreich and Lazutkin 2001]. 
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Figure 1: Left: a spiralling orbit of the explicit Euler method with the time-step e = 0.01 
applied to the Lotka-Volterra system. Right: three orbits of the Kahan's discretization 
with e = 0.1. 

3 Hirota-Kimura bases and integrals 

In this section a general formulation of a remarkable mechanism will be given, which 
seems to be responsible for the integrability of the Hirota-Kimura type (or Kahan type) 
discretizations of algebraically completely integrable systems. This mechanism is so far 
not well understood, in fact at the moment we do not know what mathematical structures 
make it actually work. 

Throughout this section / : M. n — > R n is a birational map, while hi, ifi : lR n — > R stand 
for rational, usually polynomial functions on the phase space. We start with recalling a 
well known definition. 

Definition 1 A function h : W 1 — > R is called an integral, or a conserved quantity, 

of the map f , if for every xq € R n there holds 

h(f(x)) = h(x), 

so that 

hof i (x) = h(x) VieZ. 

Convention. In the last formula and everywhere in the sequel, we use the expression 
h o f l {x) for the evaluation of the function ho f % at the point x. This is equivalent to 
h(f l (x)) and is used to spare some parentheses. 

Thus, each orbit of the map / lies on a certain level set of its integral h. As a conse- 
quence, if one knows d functionally independent integrals hi, . . . , ha of /, one can claim 
that each orbit of / is confined to an (n — cf)-dimensional invariant set, which is a common 
level set of the functions hi, . . . , h d . 
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Definition 2 A set of functions $ = (ipi, . . . ,(pi), linearly independent over M., is called 
a Hirota-Kimura basis (HK-basis) , if for every xq G W 1 there exists a vector c = 
(ci, . . . , q) ^ such that 

(c m + ... + c m ) o f(x) = Vi G Z. (8) 

For a gwen a; G M n , the set of all vectors c G M. 1 with this property will be denoted by 
K^(x) and called the null-space of the basis $ (at the point x). This set clearly is a vector 
space. 

Thus, for a HK-basis $ and for c G K$>(x) the function h = c\^>\ + ... + c\(pi vanishes 
along the /-orbit of x. Let us stress that we cannot claim that h = C\<p\ + ... + ciipi is 
an integral of motion, since vectors c G K${x) do not have to belong to K<^(y) for initial 
points y not lying on the orbit of x. However, for any x the orbit {f l (x)} is confined to 
the common zero level set of d functions 

hj = ci j Vi + • • • + cpVj = 0, j = l,...,d, 

where the vectors = (c^, . . . , cf^) G R z form a basis of K^(x). Thus, knowledge 
of a HK-basis with the null-space of dimension d leads to a similar conclusion as knowl- 
edge of d independent integrals of /, namely to the conclusion that the orbits lie on 
(n — (i)-dimensional invariant sets. Note, however, that a HK-basis gives no immediate 
information on how these invariant sets foliate the phase space W 1 , since the vectors 
and therefore the functions hj, change from one initial point x to another. 

Although the notions of integrals and of HK-bases cannot be immediately translated 
into one another, they turn out to be closely related. 

The simplest situation for a HK-basis corresponds to I = 2, dimi^<j,(x) = d = 1. In 
this case we immediately see that h = (fi/(f2 is an integral of motion of the map /. 
Conversely, for any rational integral of motion h — ipx/ipi its numerator and denominator 
(fx, if2 satisfy 

(cxcpt + c 2 (p 2 ) o f{x) = 0, % G Z, 

with c\ = 1, Ci = —h(x), and thus build a HK-basis with 1 = 2. Thus, the notion of a 
HK-basis generalizes (for I > 3) the notion of integrals of motion. 

On the other hand, knowing a HK-basis $ with dimK$(x) = d > 1 allows one to find 
integrals of motion for the map /. Indeed, from Definition [2] there follows immediately: 

Proposition 3 If $ is a HK-basis for a map f , then 

K*(f(x)) = K„(x). 

Thus, the rf-dimensional null-space K$(x) G Gr(d, I), regarded as a function of the initial 
point x G M. n , is constant along trajectories of the map /, i.e., it is a Gr(d, Z)-valued 
integral. One can extract from this fact a number of scalar integrals. 

Corollary 4 Let $ be a HK-basis for f with dim Kq>(x) = d for all x G lR n . Take a basis 
of K$(x) consisting of d vectors cW e M! and put them into the columns of a I x d matrix 
C(x). For any d-index a = (oil, . . . , a^) C {1,2, ... ,n} let C a = C ai ... ad denote the dx d 
minor of the matrix C built from the rows a\, . . . , a^- Then for any two d-indices a, f3 
the function Ca/Cp is an integral of f . 
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Proof. The functions C a are nothing other than the Grassmann-Pliicker coordinates 
of the (i-space K$(x) in the Grassmannian Gr(d,l), which are defined up to a common 
factor. More detailed, any basis of K$(f(x)) is obtained from the given basis of K$(x) via 
a right multiplication of C by a non-degenerate dxd matrix D. This yields a simultaneous 
multiplication of all C a by the common factor det D. This operation does not change the 
quotients C a /Cp. □ 
Especially simple is the situation when the null-space of a HK-basis has dimension 
d = l. 

Corollary 5 Let $ be a HK-basis for f with dimK$(x) = 1 for all x G M. n . Let K$(x) = 
[ci(x) : . . . : ci(x)] G MP' -1 . Then the functions Cj/c^ are integrals of motion for f . 

An interesting (and difficult) question is about the number of functionally independent 
integrals obtained from a given HK-basis according to Corollaries [4] and [5] We will see 
later that it is possible for a HK-basis with a one-dimensional null-space to produce more 



than one independent integral (see Theorem 13) 



The first examples of this mechanism (with d — 1) were found in [Kimura and Hirota 
2000] and (somewhat implicitly) in [Hirota and Kimura 2000]. 



4 Finding Hirota-Kimura bases 

At present, we cannot give any theoretical sufficient conditions for existence of a Hirota- 
Kimura basis $ for a given map /, and the only way to find such a basis remains the 
experimental one. Definition [2] requires to verify condition ^ for all i £ Z, which is, of 
course, impractical. We now show that it is enough to check this condition for a finite 
number of iterates f l . 

For a given set of functions $ = ((pi, . . . , cpi) and for any interval [j, k] C Z we denote 

/ ^(P(x)) .. Vl (P(x)) \ 



x \j,k]( x ) 



(9) 



V ^{f k {x)) .. <pi(f k (x)) J 
In particular, A(„ 00 00 )(x) will denote the double infinite matrix of the type (J9|. Obviously, 

ker A(_ 00)00 )(x) = K$(x). 



Thus, Definition [2] requires that dim ker -X'(_ CIO)ao ) (x) > 1. Our algorithm for detecting this 
situation is based on the following observation. 

Theorem 6 Let 

dimker X {0>s ^(x) = { 1 ~ S ^ ^ i^+t 

hold with some d for all x G M n . Then for any x G M. n there holds: 

ker A(_ 00i0o) (x) = kerX [0 ,z_d-i](a;) ) 

and, in particular, 

dimkerX^^oo)^) = d. 
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Proof. By definition, Xy jk ](x) = X\ 0)k _fl(p(x)). Therefore, applying condition (10) to 
iterates P(x) instead of x itself, we see that the kernel of any sub matrix of X(_ 0000 )(a;) 
with l — d rows, as well as the kernel of any submatrix with l — d + 1 rows, is ci-dimensional: 

dimkerX[jj +J _ d _i](x) = dim ker X^+i-^ (x) = dimker X [j+l . j+ i„ d] (x). 

Since, obviously, 

kerX[ iii+ /_ d _i](x) D ker X^j+i^ (x) C kerXy + i ii+ /_ d ](x), 

we find that all three kernels coincide, in particular, 

ker Xtfj+i^^x) = ker X [j+1J+l ^ d] (x). 

By induction, all ker Xyj+i-d-i^x), j G Z, coincide, and therefore they coincide with 
ker X(_ 00i00 )(x), as well. □ 
These results lead us to formulate the following numerical algorithm for the estimation 
of dimi£$(x) for a hypothetic HK-basis $ = (<fi, . . . , tpi). 

(N) For several randomly chosen initial points x G lR n , compute dimker X[ 0jS _i](x) for 



1 < s < I. If for every x condition (10) is satisfied with one and the same d > 1, 



then $ is likely to be a HK-basis for /, with dim K$(x) = d. 

We stress once again that generally (for general maps / and general monomial sets $) 
one will find that the I x I matrix X[o,i-i]{x) is non-degenerate for a typical x, so that 
dim K$(x) = 0. Finding (a candidate for) a HK-basis $ is a highly non-trivial task. 

Having found a HK-basis $ with dimK$(x) = d numerically, one faces the next 
problem: to prove this fact, that is, to prove that the system of equations ^ with 
i = io, z'o + 1, . . • ,io + I — d admits (for some, and then for all z £ Z) a rf-dimensional 
space of solutions. For the sake of clarity, we restrict our following discussion to the most 
important case d = 1. Thus, one has to prove that the homogeneous system 

(dip! + . . . + cupi) o f l (x) = 0, i = i ,i + l,...,i + l-l (11) 

admits for every iGM"a one-dimensional vector space of non-trivial solutions. The main 



obstruction for a symbolic solution of the system (11) is the growing complexity of the 
iterates f l {x). While the expression for f{x) is typically of a moderate size, already the 
second iterate f 2 (x) becomes typically prohibitively big. In such a situation a symbolic 



solution of the linear system (11) should be considered as impossible, as soon as f (x) is 



involved, for instance, if I > 3 and one considers the linear system with i = 0, 1, . . . , I — 1. 



Therefore it becomes crucial to reduce the number of iterates involved in (11) as far 
as possible. A reduction of this number by 1 becomes in many cases crucial! One can 
imagine several ways to accomplish this. 

(A) Take into account that, because of the reversibility f~ 1 (x,e) = f(x,—e), the neg- 
ative iterates f~ l are of the same complexity as /*. Therefore, one can reduce the 



complexity of the functions involved in (11) by choosing i = —[1/2] instead of the 
naive choice io = 0. 
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For instance, in the case I = 3 one should consider the system (11 ) with i = —1, 0, 1, and 



not with i — 0, 1, 2. However, already in the case I = 4 this simple recipe does not allow 
us to avoid considering f 2 . In this case, the following way of dealing with the system ( 11 ) 
becomes useful. 



(B) Set q = — 1 and consider instead of the homogeneous system ( 11 ) of Z equations the 
non-homogeneous system 

(c! </?i + ... + of (x) =(piof(x), i = i ,i + l,...,i + l-2, (12) 

of / — 1 equations. Having found the (unique) solution (ci(x), . . . , prove 
that these functions are integrals of motion, that is, 



c i{f{ x )) = ci(z), 



Q_i(/(x)) = C;_i(x). 



(13) 



Thus, for instance, in the case I = 4 one has to deal with the non-homogeneous system 



of equations (12) with i = —1, 0, 1. Unfortunately, even if one is able to solve this system 



symbolically, the task of a symbolic verification of eq. (13) might become very hard due 
to complexity of the solutions (ci(x), . . . , cj_i(x)J. 

This is the way taken, for instance, in [Kimura and Hirota 2000]. In that paper, the 



task of verifying the equations of the type (13) for the discrete time Lagrange top is 
performed with the following method. 

(G) In order to verify that a rational function c(x) = p(x)/q(x) is an integral of motion 
of the map x = f(x) coming from a system ([2]): 

i) find a Grobner basis G of the ideal / generated by the components of eq. (|2]), 
considered as multi-linear polynomials of 2n variables x, x of total degree 2; 

ii) check, via polynomial division through elements of G, whether the polynomial 
5(x, x) = p(x)q(x) — p(x)q(x) belongs to the ideal i". 

An advantage of this method is that neither of its two steps needs the complicated explicit 
expressions for the map /. Nevertheless, both steps might be very demanding, especially 
the second step in case of a complicated integral c(x). 



Sometimes, the task of verifying equations (13) can be circumvented by means of the 
following tricks. 



(C) Solve system (12) for two different but overlapping ranges i G [io,ia + 1 — 2] and 
i £ [ii, i\ + I — 2]. If the solutions coincide, then eq. (13) holds automatically. 



Indeed, in this situation the functions (ci(x), . . . , q_i(x)) solve the system with i e 
[zo, io + I — 2} U [ii, i\ + I — 2] consisting of more than / — 1 equations. 

A clever modification of this idea, which allows one to avoid solving the second system, 
is as follows. 



(D) Suppose that the index range i G [io,io + / — 2] in eq. (12) contains but is non- 
symmetric. If the solution of this system (ci(x, e), . . . , q_i(x, e)) is even with respect 
to e, then eqs. (13) hold automatically. 
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Indeed, the reversibility of the map / e) = f(x, — e) yields in this case that equations 



of the system (12) are satisfied for i e [— (i + I — 2),— i ], as well, and the intervals 
[io, z'o + I — 2] and [— («o + ^ — 2), — i ] overlap but do not coincide, by condition. 

Finally, the most powerful method of reducing the number of iterations to be considered 
is as follows. 

(E) Often, the solutions (ci(x), . . . , q_i(a;)) satisfy some linear relations with constant 
coefficients. Find (observe) such relations numerically. Each such (still hypothetic) 



relation can be used to replace one equation in the system (12). Solve the resulting 



system symbolically, and proceed as in recipes (C) or (D) in order to verify eqs. 



(13). 



In some (rare) cases the integrals found by this approach are nice and simple enough 



to enable one to verify eqs. (13) directly. Of course, it would be highly desirable to 
find some structures, like Lax representation, bi-Hamiltonian structure, etc., which would 
allow one to check the conservation of integrals in a more clever way, but up to now no 
such structures have been found for any of the Hirota-Kimura-type discretizations. 



5 Hirota-Kimura discretization of the Euler top 

We now illustrate the Hirota-Kimura mechanism by its application to the Euler top. This 
three-dimensional system is simple enough to enable one to perform all necessary com- 
putations symbolically, even by hand. At the same time, it provides a perfect illustration 
for many of the issues mentioned in the previous section. 

5.1 Euler top 

The differential equations of motion of the Euler top read 

xi = aix 2 x 3 , x 2 = a 2 x 3 xi, x 3 = a 3 xix 2 , (14) 

with a, being real parameters of the system. This is one of the most famous integrable 
systems of the classical mechanics, with a big literature devoted to it. We mention only 
that this system can be explicitly integrated in terms of elliptic functions, and admits 
two functionally independent integrals of motion. Actually, a quadratic function H(x) = 



~fix\ + 7 2 x| + 7 3 x| is an integral for eqs. (14), if (7, a) = 71 ai + 7 2 a 2 + 72^2 = 0. In 



particular, the following three functions are integrals of motion: 



Hi = CH3X2 — a 2 x 3 , H 2 = aix 3 — a^x^ H 3 = a 2 x 1 — a±x 2 . 
Clearly, only two of them are functionally independent because of aiHi+a 2 H 2 +a^H 3 = 0. 
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5.2 Discrete equations of motion 

The Hirota-Kimura discretization of the Euler top introduced in [Hirota and Kimura 2000] 

reads as 

xi—xi = ea 1 (x 2 X3 + x 2 x 3 ), 

x 2 - x 2 = eo 2 (x 3 xi + x 3 xi), (15) 
X3- x 3 = ea 3 (x l x 2 + xix 2 ). 

Thus, the map / : x h-> x obtained by solving ((IB"]) for x, is given by: 



1 — ect\X3 — ea±x 2 \ 
x = f(x, e) = A~ 1 (x, e)x, A(x,e) = | —£0:2X3 1 — eo 2 xi . (16) 

-603X2 —£03X1 1 J 

It might be instructive to have a look at the explicit formulas for this map: 

x\ + 2ea\x 2 X3 + e 2 xi(— 02O3X 2 + 0301X2 + a\a 2 x\) 



Xi 



Afar, el 



_ x 2 + 2eo 2 X3Xi + e 2 x 2 (a 2 a3x\ - 0301X0 + oio 2 Xo) 

X2 = w^) ' (17) 

_ X3 + 2eo3Xix 2 + e 2 X3(o 2 03xf + 03O1X2 — oio 2 x|) 
X3= A(x7e) ' 

where 

A(x, e) = det A(x, e) = 1 — e 2 {a 2 a3x\ + 0301X2 + 0102X3) — 2e 3 oio 2 03Xix 2 X3. 

As always the case for HK-type discretizations, this map is birational, and there holds 
the reversibility property: 



Apart from the Lax representation which is still missing, the discretization ( 16 ) exhibits 
all the usual features of an integrable map: an invariant volume form, a bi-Hamiltonian 
structure (that is, two compatible invariant Poisson structures), two functionally inde- 
pendent conserved quantities in involution, and solutions in terms of elliptic functions. 
The difference of its qualitative behavior as compared with non-integrable discretizations 



is striking, cf. Fig. |5.2 For further details about the properties of this discretization we 



refer to [Hirota and Kimura 2000] and [Petrera and Suris 2008]. The integrals have been 
first found in [Hirota and Kimura 2000], apparently with the help of the approach dis- 
cussed in the present work. However, since the resulting integrals are sufficiently simple 
and nice, their conservation can be easily verified by hands, therefore the paper [Hirota 
and Kimura 2000] presents them in an ad hoc form, without explaining how they have 
been discovered. We now try to reconstruct the way the results of [Hirota and Kimura 



2000] were originally found. For this aim, we apply to the map (16) the method described 
in section [3j 
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Figure 2: Left: a spiralling orbit of the explicit Euler method with the time-step e = 0.3 
applied to the Euler top. Right: a single orbit of the Hirota-Kimura discretization with the 
same time-step, lying on an invariant spatial elliptic curve (intersection of two quadrics). 



5.3 Hirota-Kimura bases 

Since all integrals of the Euler top are linear combinations of the functions x\, it is natural 
to try the set 

$ = (x?,44i) (is) 

as a HK-basis for the discrete time Euler top. An application of the numerical algorithm 
(N) suggests that the following statement holds: 



Theorem 7 The set (18) is a HK-basis for the map (16) with dimi£$(x) = 2. Therefore, 
any orbit of this map lies on the intersection of two quadrics in IR 3 . 

We will prove this theorem by finding two smaller HK-bases with d = 1 . Namely, applica- 
tion of the numerical algorithm (N) suggests that omitting any one of the four functions 
1, x\ from the basis $ leads to a HK-basis with d = 1. In other words, for every iel 3 
there exists a one-dimensional space of vectors (01,02,03) such that 

(c\x\ + c 2 x\ + C3X3) o f\x) =0, i eZ, 

as well as a one-dimensional space of vectors (d\, d 2 , c/4) such that 

{d l x\ + d 2 x\ + ck) o f(x) =0, i £ Z. 

These numerical results can be now proven analytically. 

Proposition 8 The set 

^0 — (#1) x \i ^3) 



is a HK-basis for the map (16) with dimK^ (x) = 1. At each point i6i 3 there holds: 
K^ (x) = [ci : c 2 ■ C3J = IQ>3X 2 — a 2 x 3 : ctix 3 — a^x 1 : a 2 x ± — aix 2 \- 
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Setting C3 = —1, the functions 



cAx) 



a^x\ 



a 2 x\ 



a\x\ — a 2 x\ ' 



c 2 {x) 



a\x\ — a^x\ 
oi\x\ — a 2 x\ 



(19) 



are integrals of motion of the map (16) 



Proof. We proceed according to the recipe (B), set C3 
system 



•1, and solve symbolically the 



{c x x\ + c 2 x\) o 



X, O 



0,1, 



(20) 



which involves two non-homogeneous equations for two unknowns. System (20) can be 
written as 



Cix\ + c 2 x\ 
c{x{ + c 2 x\ 



x 



3! 



(21) 



X 



3) 



where, of course, explicit formulas (17) have to be used for x^. The solution of this system 



is given by formulas (19). The components of the solution do not depend on e, therefore, 



according to the recipe (D), we conclude that functions (19) are integrals of motion of the 
map (16). □ 



It should be mentioned that the independence of the solution (ci,c 2 ) on e, or, more 
generally, the dependence through even powers of e only, which will be mentioned on 
many occasions below, starting with Proposition [9j is not granted by any well- understood 
mechanism. Rather, it is just an instance of very remarkable and miraculous cancellations 
of non-even polynomials. We illustrate this phenomenon by providing additional details 



to the previous proof. The solution of eqs. (21) by the Cramer's rule is given by ratios of 
determinants of the type 



,y>2 ,yi 2 



X 



X ■ 



Ae{a.jx1 - a i x'j)(xi + ea 1 x 2 x 3 )(x 2 + ea 2 x 3 xi)(x3 + ea^xix 2 ) 

A 2 (x, e) 



(22) 



In the ratios of such determinants everything cancels out, except for the factors 

The cancellation of the denominators A 2 (x, e) is, of course, no wonder, but the cancellation 

of the non-even factors in the numerators is rather miraculous. 

One more typical phenomenon occurs in Proposition IHJ although we have found ap- 



parently two integrals of motion (19), they turn out to be functionally dependent. Indeed, 
there holds an identity 

a x cx{x) + a 2 c 2 (x) = a 3 , 

so that for each 1GR 3 the space K$ (x) is orthogonal to the constant vector (a 1; a 2 , a 3 ). 
If one would have guessed this relation numerically, one could simplify the computation 
of the integrals C\, c 2 by considering the system 



c\x\ + c 2 x\ 
C1OL1 + c 2 a 2 



X 



3) 
"3, 



(23) 



instead of (21 ). Observe that existence of a linear relation allows one to reduce a number 



of iterates of / involved in the linear system (in the present situation, the system (23) 
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contains no iterates of / at all!). The latter system would lead to the same formulas (19) 



however, in this case one could not argue as in (D), and would be forced to prove that 



the functions (19) are integrals of motion directly, by verifying for them equations (13). 



Anyway, the existence of the HK-basis $o yields existence of only one independent 
integral of the map /, which is not enough to assure the integr ability of /. 

Proposition 9 The set 

$i = {xl,x 2 2 ,l) 

is a HK-basis for the map (16) with dim K^ 1 (x) = 1. At each point i6i 3 there holds: 

K^(x) = [di : d 2 : — 1], 

where 



di(x) 



0:2(1 — €0301X2) 
a 2 x\ — a\x\ 



d 2 (x) 



Oi(l — e a 2 a 3 x\) 
Oix| — a 2 xf 



(24) 



These functions are integrals of motion of the map (16) 



Proof. Following again prescription (B), we set d± 
non-homogeneous system 

{dix\ + d 2 x\) o f(x) = 1, 



1, and solve symbolically the 



0,1, 



or 



d\x\ + d 2 x\ 
d{x\ + d 2 x\ 



1, 
1. 



The solution is given by eq. (24), due to eq. (22) and 



X; 



4eOj(l - e 2 a j a k x 2 i )(xi + eaix 2 x^)(x 2 + ea 2 x 3 x 1 )(x 3 + ea 3 X!X 2 ) 

A 2 (x,e) 



This time its components do depend on e, but are manifestly even functions of e. Every- 
thing non-even luckily cancels, again. Therefore, the argument (D) is still applicable, so 
that the functions (24) are integrals of motion of the map /. □ 



Functions (24) are again functionally dependent, because of 

OLld\(x) + 02^2(2) = e 2 0i0203. 

However, they are, clearly, functionally independent on the previously found functions 



(19), because ci,c 2 depend on x 3 , while di,d 2 do not. 

Of course, the permutational symmetry yields that each of the sets of monomials $2 



{x 2 , 1) and $3 = (x 2 , x§, 1) is a HK-basis, as well, with dimif$ 2 



. x 



dim K§ 3 (x) 



Any two of the four found one-dimensional null-spaces span the full null-space K<&{x). In 
particular, Kq, (x) lies in K^^x) © K^ 2 (x). 

Summarizing, we have found a HK-basis with a two-dimensional null-space, as well as 
two functionally independent conserved quantities for the HK-discretization of the Euler 
top. Both results yield integrability of this discretization, in the sense that its orbits 
are confined to closed curves in M 3 . Moreover, each such curve is an intersection of two 
quadrics, which in the general position case is an elliptic curve. 



14 



6 Hirota-Kimura-type discretization of the Clebsch 
system 

6.1 Clebsch system 

The motion of a rigid body in an ideal fluid can be described by the so called Kirchhoff 
equations [Kirchhoff 1870]: 

OH OH 

m — m x — hp x — — , 

dm dp (25) 

dH 

P = P X ^— > 
am 

with H being a quadratic form in m = (mi, m 2 ,m 3 ) G M 3 and p = (pi,P2,P3) £ 
here x denotes vector product in M 3 . The physical meaning of m is the total angular 



momentum, whereas p represents the total linear momentum of the system. System (25) 
is Hamiltonian with the Hamilton function H(m,p), with respect to the Poisson bracket 

{m^rrij} = m k , {m^pj} = p k , (26) 

where (i,j,k) is a cyclic permutation of (1,2,3) (all other pairwise Poisson brackets of 
the coordinate functions are obtained from these by the skew-symmetry, or otherwise 
vanish). A detailed introduction to the general context of rigid body dynamics and its 
mathematical foundations can be found in [Marsden and Ratiu 1999]. 

A famous integrable case of the Kirchhoff equations was discovered in [Clebsch 1870] 
and is characterized by the Hamilton function H — | Y^=i( m i '+ u iPi) ■ The corresponding 
equations of motion read: 

rh = p x Qp, 
p = p x m, 

where Q = diag(a>i, o>2, 0*3) is the matrix of parameters, or in components: 

mi = (UJ 3 - UJ 2 )P2P3, 

m 2 = (uji - co 3 )p 3 pi, 
m 3 = (uj 2 -Ui)piP2, 

pi = m 3 p 2 - m 2 p 3 , 

p 2 = m x p 3 - m 3 pi, 

P3 = -m 2 p x - m x p 2 . 

This is the system which will be called the Clebsch system hereafter. For an embedding of 
this system into the modern theory of integrable systems see [Perelomov 1990,Reyman and 
Semenov-Tian-Shansky 1994]. The Clebsch system possesses four independent quadratic 
integrals: 

Hi = m\ + m\ + mj + uipl + uj 2 p\ + uJ 3 p\, (27) 
H 2 = ujim\ + uo 2 m\ + uo 3 m\ - uo 2 uo 3 p\ - 0J 3 u x p\ - LOiUJ 2 p\, (28) 

#3 = Pl+P 2 2+Pl (29) 
H 4 = mipi + m 2 p 2 + m 3 p 3 . (30) 
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These integrals are in involution with respect to the bracket (26), moreover, H 3 ,H^ are 
its Casimir functions (are in involution with any function on the phase space). However, 
the Hamiltonian structure will not play any role in the present paper. The set of linear 
combinations of the quadratic Hamiltonians Hi, if 2, H 3 coincides with the set of linear 
combinations of the functions 



h = vl + 
h = Pl + 



9 9 

mi mi 





+ 




U) X - u 3 






2 








+ 


m\ 






lo 2 - UJ 3 




LU 2 — UJi 


2 






mi 


+ 


m 2 2 



h = Pt + 

uj 3 —uj 2 uj 3 - Vi 



For instance 



Hi = uili + uj 2 h + CO3I3, Hi = —CU2U3I1 — U3UJ1I2 — ujiuj2h, H 3 = I\ + J 2 + h- 



6.2 Discrete equations of motion 

Applying the Hirota-Kimura (or Kahan) approach to the Clebsch system, we arrive at 
the following discretization, proposed in [Ratiu 2006]: 

rhi - mi = e(^3 - ^2X^3 + P2P3), 

m 2 - m 2 = e(u)i - u 3 )(p 3 pi + p 3 pi), 

rh 3 -m 3 = t{u 2 - u)i){p~ip 2 + P1P2), 

Pi- Pi = t{rh 3 p 2 + m 3 p 2 ) - e(m 2 p 3 + m 2 p 3 ), 

P2-P2 = t(rhip 3 + mip 3 ) - e(fh 3 pi + m 3 pi), 

P3-P3 = e(m 2 pi + m 2 pi) - t{mip 2 + mip 2 ). 

In matrix form this can be put as 

g) = Q. 

where 

M(m,p, e) 



and the abbreviation = 
the b national map / : IR 6 

(^=f(m,p,e) = M-\m, P ,e)(^y (31) 



/ 1 











tu 23 p 3 


€UJ23P2\ 





1 





e^3iP3 





eu 3 ipi 








1 


e^i2P2 


eui 2 pi 








eP3 




1 


-em 3 


em 2 







epi 


em 3 


1 


—emi 


V eP2 


-e P i 





-em 2 


emi 


1 / 



uj,i — ujj is used. The solution of this 6x6 linear system yields 
-> M 6 , 
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called hereafter the discrete Clebsch system. As usual, the reversibility property holds: 

f' 1 (m,p,e) = f(m,p,-e). (32) 

A remark on the complexity of the iterates of / is in order here. Each component of 
(fh,p) = f(m,p) is a rational function with the numerator and the denominator being 
polynomials on m^^pk of total degree 6. The numerators of pk consist of 31 monomials, 
the numerators of rhk consist of 41 monomials, the common denominator consists of 28 
monomials. It should be taken into account that the coefficients of all these polynomials 
depend, in turn, polynomially on e and u>k, which additionally increases their complex- 
ity for a symbolic manipulator. Expressions for the second iterate swell to astronomical 
length prohibiting naive attempts to compute them symbolically. Using MAPLE's Large- 
Expressions package [Carette et al. 2006] and an appropriate veiling strategy it is however 
possible to obtain f 2 (m,p) with a reasonable amount of memory. Some impression on 
the complexity can be obtained from Table [1} The resulting expressions are too big to 
be used in further symbolic computations. Consider, for instance, the numerator of the 
Pi-component of f 2 (m,p). As a polynomial of rrik,Pk, it contains 64 056 monomials; their 
coefficients are, in turn, polynomials of e and u^, and, considered as a polynomial of the 
phase variables and the parameters, this expression contains 1 647 595 terms. 





deg 


de s P1 


de s P2 


de § P3 


deg mi 


deg m2 


deg m3 


Common denominator of f 2 


27 


24 


24 


24 


12 


12 


12 


Numerator of pi-comp. of f 2 


27 


25 


24 


24 


12 


12 


12 


Numerator of p2-comp. of f 2 


27 


24 


25 


24 


12 


12 


12 


Numerator of p3-comp. of f 2 


27 


24 


24 


25 


12 


12 


12 


Numerator of m r comp. of f 2 


33 


28 


28 


28 


15 


14 


14 


Numerator of m2-comp. of f 2 


33 


28 


28 


28 


14 


15 


14 


Numerator of m3-comp. of f 2 


33 


28 


28 


28 


14 


14 


15 



Table 1: Degrees of the numerators and the denominator of the second iterate f 2 (m,p) 



6.3 Phase portrait and integrability 

We now address the problem whether the discrete Clebsch system is integrable. Figs. 



|3| and |4j show plots of the discrete Clebsch system (31), produced with MATLAB, for 
two different sets of parameters values. These plots indicate a quite regular behavior of 
the orbits of the discrete Clebsch system. Each orbit seems to fill out a two-dimensional 
surface in the 6-dimensional phase space. Leaving aside the Hamiltonian aspects of in- 



tegrability, we are interested just in this simpler issue: do orbits of the map (31) lie 
on two-dimensional surfaces in IR 6 ? A usual way to establish such a property would be 
to establish the existence of four functionally independent conserved quantities for this 
map. (We note in passing that plots of orbits are not very reliable in deciding about 
integrability. For instance, there are indications that the Kahan's discretization ^ of the 
Lotka-Volterra system is non-integrable, even if its orbits visually lie on closed curves in 
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(a) m 1 ,m 2 ,m 3 (b) p 1 ,p 2 ,p 3 

Figure 3: An orbit of the discrete Clebsch system with lo\ = 1, uj 2 = 0.2, uj 3 = 30 and 
e = 1; initial point (m ,Po) = (1, 1, 1, 1, 1, 1). 

the phase plane. A strong magnification unveils the existence of very small regions in the 
phase plane with a chaotic behavior.) 

We will show that the answer to the above question is in affirmative. For this aim, 
we apply the approach based on the notion of HK-basis. As a first step, we apply the 
numerical algorithm (N) to the maximal set of monomials, which includes all monomials 



of which the integrals (27)-(30) of the continuous Clebsch system are built: 



<pi(m,p) = pi, tp2(m,p) = p 2 , ^{m,p) = P^ 

V9 4 (m,p) = m\, <P5(m,p) = m\, <^ 6 (m,p) = m\, 

(p 7 (m,p) = mipi, (f 8 (m,p) = m 2 p 2 , (pg(m,p) = m 3 p 3 , 

(pio{m,p) = 1. 

We come to the following result: 
Theorem 10 The set of functions 

$ = (piiPiiPlirnl,ml,ml,m 1 p 1 ,m 2 p 2 ,m 3 p 3 , 1) 



is a HK-basis for the map (31), with dimi^,j,(m,p) = 4. Thus, any orbit of the map (31) 
lies on an intersection of four quadrics in IR 6 . 



At this point Theorem 10 remains a numerical result, based on the algorithm (N). A 
direct symbolical proof of this statement is impossible, since it requires dealing with f\ 
i G [—4,4], and the fourth iterate f A is a forbiddingly large expression. In order to prove 
Theorem [lO] and to extract from it four independent integrals of motion, it is desirable to 
find HK-(sub)bases with a smaller number of monomials, corresponding to some (prefer- 
ably one-dimensional) subspaces of K$(m,p). A much more detailed information on the 
HK-bases is provided by the following statement. 
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(a) mi,m 2 ,m 3 



(b) Pl,P2,P3 



Figure 4: An orbit of the discrete Clebsch system with u\ = 0.1, uo 2 = 0.2, uj 3 = 0.3 and 
e = 1; initial point (mo,Po) = (1, 1, 1, 1, 1, !)■ 



Theorem 11 The following four sets of functions are HK-bases for the map (31) with 
one- dimensional null-spaces: 



$ 

$ 3 



(pl,pl,pl,ml,ml,ml : m 2 p 2 ) , 
(Pi,pI,pI, m l, m l, m h m 3 p 3 ). 



(33) 
(34) 
(35) 
(36) 



If all the null-spaces are considered as subspaces of 1R ; so that 



K ®0 


= [ct: 


c 2 : 


c 3 : 


: 


: 


: 


: 


: 


: 


Cio], 




= [aii ■ 


ct 2 ■ 


a 3 : 


«4 : 


a 5 : 


«6 : 


«7 : 


: 


: 


o], 


K$ 2 


= \Pi- 


(32: 


A = 


ft: 




A>: 


: 




: 


0]. 


K *3 


= hi ■ 


72 : 


7s : 


74 : 


7s : 


76 : 


: 


: 


79 : 


0]. 



then there holds: 



K, 



K, 



<p 2 



Also this statement was first found with the help of numerical experiments based on the 
algorithm (N). In what follows, we will discuss how these claims can be given a rigorous 
(computer assisted) proof, and how much additional information (for instance, about 



conserved quantities for the map (31 )) can be extracted from such a proof. 
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6.4 First HK-basis 



Theorem 12 The set (33) is a HK-basis for the map (31) with dimfT$ (m,p) = 1. 
each point (m,p) G M 6 there holds: 

K^ (m,p) = [ci : c 2 : c 3 : c w ] 

1 + e 2 (uj 1 - lu 2 )pI + e 2 (uj 1 - cu 3 )pl ^ 1 + e 2 (tu 2 ~ ^i)p\ + e 2 (^2 - ^3)^3 . 



At 



pI + p 2 2 +pI 

1 + e 2 (u 3 - u x )p\ + e 2 (u 3 - u 2 )p 2 

vl+vl+vl 

— + e 2 Ui : - + e 2 uj 2 '■ j + e 2 ^3 : -1 



Pi + v\ + pi 



where 



J(m,p,e) 



P1+P2+ PI 



1 - e 2 (u;ip? + io 2 p\ + u 3 pl) ' 



The function (38) is an integral of motion of the map {Sly . 

Proof. The statement of the theorem means that for every (m,p) G 
solutions of the homogeneous system 



(37) 



(38) 



the space of 



(cip 2 + c 2 p\ + C3P3 + cio) o f(m,p) = 0, 



0. 



is one-dimensional. This system involves the third iterate of /, therefore its symbolical 
treatment is impossible. According to the strategy (B), we set cio = — 1 and consider the 
non-homogeneous system 



{cip\ + c 2 p\ + C3P3) o f(m,p) = 1, 



0,1,2. 



(39) 



This system involves the second iterate of /, which still precludes its symbolical treatment. 
There are now several possibilities to proceed. 

• First, we could follow the recipe (E) and find further information about the solutions 
Ci. For this aim, we plot the points {ci{m ) p) 1 c 2 {m 1 p),c 3 {m 1 p)) for different initial 
data (m,p) G M 6 . Fi gure pi shows such a plot, with 300 initial data {m,p) randomly 
chosen from the set [0, l] 6 . The points (ci(m,p), c 2 (m,p), c 3 (m,p)) seem to lie on a 
line in IR 3 , which means that there should be two linear dependencies between the 
functions Ci,c 2 and c 3 . In order to identify these linear dependencies, we run the 
PSLQ algorithm [Ferguson and Bailey 1991, Ferguson, Bailey and Arno 1999] with 
the vectors (ci,c 2 , 1) as input (see Remark after the end of the proof, concerning 
implementation of this step). On this way we obtain the conjecture 

ci - c 2 = e 2 (^i - u 2 ). 

Similarly, running the PSLQ algorithm with the vectors (c 2 , c 3 , 1) as input leads to 
the conjecture 

c 2 - c 3 = e 2 (u 2 - uj 3 ). 
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Komp. 1 ,2.3 




-12 _,, -15 



Figure 5: Plot of the coefficients ci, c 2 , C3 



Having identified (numerically!) these two linear relations, we use them instead 
of two equations in the system (39), say the equations for i = 1,2. The resulting 
system becomes extremely simple: 



cipj + c 2 pj + c 3 p 2 3 = 
c\ ~ c 2 = e 2 (^i - u 2 ) 
c 2 - c 3 = e 2 (u 2 - u 3 ) 



It contains no iterates of / at all and can be solved immediately by hands, with the 



result (37). It should be stressed that this result still remains conjectural, and one 



has to prove a posteriori that the functions c±, c 2 , C3 are integrals of motion. 

Alternatively, we can combine the above approach based on the prescription (E) 
with the recipe (D). For this, we use just one of the linear dependencies found 
above to replace the equation in (39) with % — 2, and then let MAPLE solve the 



remaining system. The computation takes 22,33 sees, on a 1.83 Ghz Core Duo 



PC and consumes 32,43 MB RAM. The output is still as in (37), but arguing this 



way one does not need to verify a posteriori that ci,c 2 ,Cs are integrals of motion, 
because they are manifestly even functions of e, while the symmetry of the linear 
system with respect to e has been broken. 

To finish the proof along the lines of the first of the possible arguments above, we show 



how to verify the statement that the function J in (38) is an integral of motion, i.e., that 



p\ + p\ + pi 



pI + pI + pI 



1 - e 2 (cjipf + u 2 p\ + u 3 pl) 1 - e 2 (wip? + u 2 pl + cj 3 p|) ' 
This is equivalent to 

~2 2 1 ~2 2 1 ~2 2 

P1-P1+P2-P2+P3- Pz 
= e 2 [{u 2 - uji)(pIp 2 2 - pipl) + (w 3 - ^ 2 )&3 - P3P2) + (wi - w 3 )(p!p? 



~2 2\ 
PlP3) 
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On the left-hand side of this equation we replace pi — pi through the expressions from the 



last three equations of motion (31 ), on the right-hand side we replace e{uj k —ujj){pjp k +pjp k ) 



by rrii — mi, according to the first three equations of motion (31 ). This brings the equation 
we want to prove into the form 

(pi + Pi)(m 3 p 2 + m 3 p 2 - m 2 p 3 - m 2 p 3 ) + 
(pa + P2)(m 1 p 3 + mip 3 - m 3 p x - m 3 pi) + 
(p~3 + Pz){m 2 pi + m 2 pi - m x p 2 - m x p 2 ) = 

= (piP2 - PiP2){rh 3 - m 3 ) + (p 2 p 3 - P2P3)(m 1 - mi) + (p 3 pi - p 3 pi)(fh 2 - m 2 ). 

But the latter equation is an algebraic identity in twelve variables m k ,p k ,rh k ,p k . This 
finishes the proof. □ 



Remark In the above proof and on many occasions below we make use of the PSLQ 
algorithm in order to identify possible linear relations among conserved quantities. Its 
applications are well documented in the literature on Experimental Mathematics [Borwein 
and Bailey 2003, Borwein, Bailey and Girgensohn 2004], so that we restrict ourselves here 
to a couple of minor remarks. We apply the PSLQ algorithm to the numerical values 
of (the candidates for) the conserved quantities obtained from the algorithm (N). We 
note that it is crucial to apply the PSLQ algorithm with many different initial data; 
from the large amount of possible linear relations one should, of course, filter out those 
relations which stay unaltered for different initial data. It proved useful to perform these 
computations with rational data (initial values of phase variables and parameters of the 
map) as well as with high precision floating point numbers. In our experiments we have 
been able to automate this task to a large extent. All computations of this kind were 
performed on an Apple MacBook with a 1.83 GHz Intel Core Duo processor and 2 GB of 
RAM. 



6.5 Remaining HK-bases 

We now consider the remaining HK-bases $i,$2 and $3. Here we are dealing with the 
three linear systems 

(«iPi + a 2 p\ + a 3 p\ + a 4 m\ + a 5 ml + a e ml) o f l (m,p) = m\P\ o f(m,p), (40) 
(Pipl + fcpl + fcpl + faml+feml + fom^of (m,p) = m 2 p 2 o f(m,p), (41) 
(7iPi + 72P2 + 73P3 + 74?™? + 75^2 + 76^1) f l (m,p) = m 3 p 3 o f l (m,p), (42) 

already made non-homogeneous by normalizing the last coefficient in each system, as in 
recipe (B), with I = 7. The claim about each of the systems is that it admits a unique 
solution for i G Z. It is enough to solve each system for two different but intersecting 
ranges of / — 1 = 6 consecutive indices i, such as i £ [— 2,3] and i e [—3,2], and to 
show that solutions coincide for both ranges (recipe (C)). Actually, since the index range 
i G [—2,3] is non-symmetric, it would be enough to consider the system for this one range 
and to show that the solutions atj, /3j, jj are even functions with respect to e (recipe (D)). 
However, symbolic manipulations with the iterates f l for i = ±2, ±3 are impossible. In 
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what follows, we will gradually extend the available information about the coefficients 
o.j , f3j , jj , which at the end will allow us to get the analytic expressions for all of them 
and to prove that they are integrals, indeed. 

6.6 First additional HK-basis 



Theorem 11 shows that, after finding the HK-basis $o with dim K$ Q (x) = 1 it is enough 
to concentrate on (sub)-bases not containing the constant function ifio(m,p) = 1. It turns 
out to be possible to find a HK-basis without <p w and with a one- dimensional null-space, 
which is more amenable to a symbolic treatment than $ 1 ,$ 2 ,^ ) 3- Numerical algorithm 
(N) suggests that the following set of functions is a HK-basis with d — 1: 

* = (pi,pi,pl,m 1 p 1 ,m 2 p 2 ,m 3 p 3 ). (43) 



Theorem 13 The set (43) is a HK-basis for the map (31) with dim K^(jn,p) = 1. At 
every point (m,p) 6 M 6 there holds: 

K-$(m,p) — [— 1 : — 1 : — 1 : d 7 : d$ : dg), 

with 

( p 2 + 2 + 2 )(1 + e 2 d (2) + £ 4 d (4) + ^(6)) 

d k = — — — — = 5_i ; k = 7, 8, 9, (44) 

A = mm + m 2 p 2 + m 3 p 3 + e 2 A (4) + e 4 A (6) + e 6 A (8) , (45) 

where and are homogeneous polynomials of degree 2q in phase variables. In 

particular, 

df' = ml + ml + ml + (u 2 + u 3 - 2ui)pl + (u 3 - u 2 )pl + (u 2 - u 3 )pl, 
<ig 2) = ml + ml + ml + (oo 3 - uoi)pl + (uj 3 + tt\ — 2uj 2 )p\ + (ui - oo 3 )pl, 
<4 2) = ml + ml + ml + (oo 2 - uJi)pl + (ui - uo 2 )pl + (u)\ + uo 2 - 2uo 3 )pl, 

and 

AW = m lPl d? + m 2V2 df + m m 4 2 \ 

(All other polynomials are too messy to be given here.) The functions dj, dg, dg are integrals 
of the map (31). They are dependent due to the linear relation 

(oo 2 - uj 3 )d 7 + (uj 3 - uJi)d 8 + (uj 1 - u 2 )dg = 0. (46) 

Any two of them are functionally independent. Moreover, any two of them together with 
J are still functionally independent. 

Proof. As already mentioned, numerical experiments suggest that for any (m,p) e IR 6 
there exists a one-dimensional space of vectors (di, d 2 , d 3 , d 7 , d$, dg) satisfying 

(dipl + d 2 p 2 2 + d 3 p 2 3 + d 7 m x pi + d s m 2 p 2 + d 9 m 3 p 3 ) o f(m,p) = 
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for i = 0,1,.. .,5. According to recipe (A), one can equally well consider this system 
for i = —2, —1, . . . , 3, which however still contains the third iterate of / and is therefore 
not manageable. Therefore, we apply recipe (E) and look for linear relations between the 
(numerical) solutions. Two such relations can be observed immediately, namely 

d 1 = d 2 = d 3 . (47) 

Accepting these (still hypothetical) relations and applying recipe (B), i.e., setting the 



common value of (47) equal to —1, we arrive at the non-homogeneous system of only 3 
linear relations 

{d 7 mipi + d 8 m 2 p 2 + d 9 m 3 p 3 ) o f (m, p) = (pj + pj + p 2 3 ) ° P (m, p) (48) 

for i = —1,0,1. Fortunately, it is possible to find one more linear relation between 
d 7 ,d 8 ,d 9 . This was discovered numerically: we produced a three-dimensional plot of the 
points {df{m,p), d%(m,p), dg(m,p)) which can be seen in Fig.[6]in two different projections. 
This figure suggests that all these points lie on a plane in IR 3 , the second picture being a 



Komp. 4,5.6 Komp. 4.5.6 




Figure 6: Plot of the points (dr,ds,dg) for 729 values of (m,p) from a six-dimensional 
grid around the point (1, 1, 1, 1, 1, 1) with a grid size of 0.01 and the parameters e = 0.1, 
oj\ = 0.1, tjj 2 = 0.2, UJ3 = 0.3. 



"side view" along a direction parallel to this plane. Thus, it is plausible that one more 
linear relation exists. With the help of the PSLQ algorithm this hypothetic relation can 
be then identified as eq. (46). Now the ansatz (48) is reduced to the following system of 
three equations for (d?, d 8 , dg), which involves only one iterate of the map /: 



{d-jmxpx + d 8 m 2 p 2 + d 9 m 3 p 3 ) o f l (m,p) = (pf+pl+pj) o /*(m,p), 

(co> 2 - ^3)^7 + (co> 3 - LUl)d S + (o> 2 - ^2)^9 = 0. 



0,1, 



(49) 



This system can be solved by MAPLE, resulting in functions given in eqs. (44), (45) 



These (long) expressions can be found in [Worksheets] . They are manifestly even functions 
of e, while the system has no symmetry with respect to e ^ — e. This proves that they 
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are integrals of motion for the map /. This argument slightly generalizes the recipes (D) 
and (E), and, since it is used not only here but also on several further occasions in this 
paper, we give here its formalization. 

Proposition 14 Consider a map f : 1R 6 — > 1R 6 depending on a parameter e, reversible in 
the sense of eq. (32). Let I(m,p, e) be an integral of f , even in e, and let Ai, A 2 , A3 e 1R. 
Suppose that the set of functions $ = (<pi, . . . , ^4) is such that the system of three linear 
equations for (ai, a 2 , a 3 ), 

(aiifi! + a 2 ip 2 + a 3 ip 3 ) o f (m,p,e) = y> 4 ° /* (m,p,e), i = 0,1, ^ 
Aiai + A 2 a 2 + A 3 a 3 = I(m,p, e), 

admits a unique solution which is even with respect to e. Then this solution (01,02,03) 
consists of integrals of the map f , and $ is a HK-basis with dim K<p(m,p) = 1. 



Proof. Since (01,02,03) are even functions of e, they satisfy also the system (50) with 
e 1— > — e, which, due to the reversibility (32), can be represented as 

(tii^i + a 2 cp 2 + a 3 (p 3 ) o f l (m,p,e) = y? 4 o f l (m,p,e), % = 0, -1, ^ 
Aia x + A 2 a 2 + A 3 a 3 = I(m,p,e). 



Since the functions (01, a 2} a 3 ) are uniquely determined by any of the systems (50) or (51 ), 
we conclude that they remain invariant under the change (m,p) 1— > f(m,p, e), or, in other 
words, that they are integrals of motion. Finally, we can conclude that these functions 
satisfy equation (01^1 + a 2 ip 2 + a 3 ip 3 ) o f % = <yj 4 o f % for all 4 G Z (and can be uniquely 
determined by this property), and that linear relation A\d\ + A 2 a 2 + A 3 a 3 = I is satisfied, 
as well. □ 

Application of Proposition 14 to system (49) shows that 07, ds, dg are integrals of mo- 
tion, since they are even in e. Note that here, as always in similar context, the evenness 
of solutions is due to "miraculous cancellation" of the equal non-even polynomials which 
factor out both in the numerators and denominators of the solutions. In the present 
computation, these common non-even factors are of degree 2 in e. 

It remains to prove that any two of the integrals dj,ds, dg together with the previously 
found integral J are functionally independent. For this aim, we show that from such 
a triple of integrals one can construct another triple of integrals which yields in the 
limit e — > three independent conserved quantities H 3 ,H 4 , Hi of the continuous Clebsch 
system. Indeed: 

j = p l +p l +p l + (e 2 )=H 3 + 0(e 2 ), 

- — = m 1 p 1 +m 2 p 2 + m 3 p 3 + 0(e 2 ) = H^ + 0(e 2 ). 
dk+e 

On the other hand, it is easy to derive: 

and, taking this into account and computing the terms of order e 4 , one finds: 
d 



L }-l-e 2 {u 2 -uj l )J = e\ 
from which one easily extracts H\. This proves our claim. □ 



1 - e 2 (iu 2 - wi) J = e\u 2 - uj x ){m{ + cu 2 H 2 3 - 2# 3 # 1 ) + 0(e b ), 

da 
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Remark With the basis we encounter for the first time the following interesting 
phenomenon: it can happen that a HK-basis with a one- dimensional null-space provides 
several (in this case two) functionally independent integrals. With Theorem 13, we estab- 
lished existence of three independent conserved quantities and two HK-bases with linearly 
independent null-spaces. So, every orbit of the discrete Clebsch system is shown to lie 
in a three-dimensional manifold which belongs to an intersection of two quadrics in IR 6 . 
The aim of the following is to find one more independent integral and two more HK-bases 
with one-dimensional null-spaces linearly independent on K$ , Kq,. 



6.7 Second additional HK-basis 



From the (still hypothetic) properties (40)-(42) of the bases $1, $2, $3 there follows that 
for any (m,p) G M 6 the system of linear equations 

{giPl+g2Pl+g3Pl+gml+g5ml+g 6 ml)of(m,p) = {mipi+m 2 p 2 +m 3 p 3 )op{m,p) (52) 

has a unique solution (g> 1; g 2 , g 3 , g±, g$, g$). Indeed, the solution should be given by 

g - = aj + (3j + 7y, j = 1, . . . , 6. (53) 



As for the bases &i,Q 2 ,Q 3 , the solution of (52) can be determined by solving these 
equations for two different but intersecting ranges of 6 consecutive values of i, say for 
i G [—3,2] and i G [—2,3]. However, it turns out that, due to the existence of several 
linear relations between the solutions gj, system (52) is much easier to deal with than 
systems (40)-(42), so that the functions gj can be determined and studied independently 
of a j ,P j ,'Yj. 

Theorem 15 The set of functions 

© = (pi,pi,pl,ml,ml,ml,mipi + m 2 p2 + rri3p 3 ) 



is a HK-basis for the map (31) with dim K®(m, p) = 1. At every point (m,p) G M 6 there 
holds: 

Ke(m,p) = [gi : g 2 ■ #3 : ■ 9s ■ 9& ■ -!]• 
Here gi,g 2 ,g 3 are integrals of the map f31\) given by 



9k 



2(p? + P2 + Ps)A 



k — 1 . 2 ; 3 ; 



where g^ are homogeneous polynomials of degree 2q in phase variables, and A is given 
in eq. (45). For instance, 

(4) 

9h' 

Integrals flag's, #6 are given by 

92 ~ #3 



1H\ — H 3 H 1 + Ld k Hn . 



9i 



oj 2 - 0J 3 



9r> 



93 - 9i 

L0 3 - LOi 



96 



9i ~ 92 

UJl - IQ 2 
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Proof. Since system (52) involves too many iterates of / for a symbolical treatment, we 



look for linear relations between the (numerical) solutions of this system. Application of 



the PSLQ algorithm allows us to identify three such relations, as given in eq. (54). This 



reduces system (52) to the following one: 



9i Pi + 



+ 

UJ\ — UJ 3 UJ\ — LU 2 



772 



+ 

^ 3 - UJ 2 0J Z — LUi 



92[P 2 + 

\ Ul 2 — ^3 UJ 2 — UJ\, 

o f(m,p) = (mipi + m 2 p 2 + m 3 p 3 ) o f(m,p). 



(54) 



Thus, one can say that we are dealing with a reduced Hirota-Kimura basis consisting of 
/ = 4 functions 

e = (h,i 2 ,i 3 ,H i ), 



see (31). Interestingly, this is a basis of integrals for the continuous-time Clebsch system, 



but we do not know whether this is just a coincidence or has some deeper meaning. System 
(54) has to be solved for two different but intersecting ranges of I — 1 = 3 consecutive 



indices i. It would be enough to show that the solution for one non-symmetric range, e.g., 
for % e [0, 2], consists of even functions of e. However, this non-symmetric system involves 
with necessity the second iterate f 2 . To avoid dealing with f 2 , one more linear relation 
for gi,g 2 ,g 3 would be needed. Such a relation has been found with the help of PSLQ 
algorithm, it does not have constant coefficients anymore but involves the previously 
found integrals d7,d$,d$: 



(u 2 - u 3 )g 1 + (u 3 - u 1 )g 2 + (u t - u 2 )g 3 = -(u 2 - ^ 3 )(^3 - Ui)(ds - d 7 ). 



(55) 



Of course, due to eq. (46), the right-hand side of eq. (55) can be equivalently put as 



^{u 3 - wi)(wi - uj 2 ){d 9 - d s ) = -(loi - io 2 ){u 2 - u 3 )(d 7 - d 9 ). 



The linear system consisting of eq. (54) for i — 0, 1 and eq. (55) can be solved by MAPLE 



with the result given in theorem. Since (dr, d%, dg) are already proven to be integrals of 



motion, and since the solutions (gi, g 2 , g 3 ) are manifestly even in e, Proposition 14 yields 
that (pi, g 2 , g 3 ) are integrals of the map /. □ 
Theorem [15] gives us the third HK-basis with a one-dimensional null-space for the 



discrete Clebsch system. Thus, it shows that every orbit lies in the intersection of three 
quadrics in M 6 . What concerns the integrals of motion, it turns out that the basis does 
not provide us with additional ones: a numerical check with gradients shows that integrals 
5'i)5 , 2,5 , 3 are functionally dependent from the previously found ones. At this point we are 
lacking one more HK-basis with a one-dimensional null-space, linearly independent from 
K$ , Ky, K e , and one more integral of motion, functionally independent from J and 
G?7, ds- 



6.8 Proof for the bases $i,$2, $3 



Now we return to the bases $i, $2; ^3 discussed in Sect. 6.5 In order to be able to 



solve systems (40)-(42) symbolically and to prove that the solutions oij,Pj,~/j are indeed 
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integrals, we have to find additional linear relations for these quantities (recipe (E)). 
Within each set of coefficients we were able to identify just one relation: 



(ui - uj 3 )a 5 
(uj 2 - u 3 )p 4 
(u 3 - o; 2 )74 



(ui - u 2 )a & , 

(u 2 - UJl)Pe, 

(u 3 - Wi)7 5 . 



(56) 
(57) 
(58) 



This reduces the number of equations in each system by one, which however does not 
resolve our problems. A way out consists in looking for linear relations among all the 
coefficients acj, Pj,jj- Remarkably, six more independent linear relations of this kind can 
be identified: 

a 4 = Ps = 7e, (59) 



(60) 
(61) 



a 2 — a 3 — (u; 2 — tj 3 )a 4 p 2 — 


Ps~ 


(a? 2 - 


OJz) Pi 


_ 72 - 


73 - 


(uj 2 - 


^3)74 


to 2 — UJ 3 

«3 - oti - (u 3 - ujx)a 5 _ P 3 - 


L0 3 

Pi - 


- Ui 

(u 3 - 




_ 7s - 


UJt 

7i - 


— <jJ 2 

{u) 3 - 


1 

^1)75 


C0 2 — UJ 3 


UJ 3 


- Ui 








— UJ 2 




are two more similar relations: 
















ai — a 2 — (uji — co 2 )a 6 Pi - 


■P2- 


(ujx - 


- uj 2 )Pe 


_ 7i - 


■ 72 - 


{u)\ - 


- a; 2 )7 6 


UJ 2 - to 3 


U 3 


- Ui 








- UJ 2 





but they follow from the already listed ones (56)-(61). We stress that all these linear 
relations were identified numerically, with the help of the PSLQ algorithm, and remain 
at this stage hypothetic. 



With nine linear relations (56)-(61), we have to solve systems (40)-(42) simultaneously 
for a range of 3 consecutive indices i. Taking this range as i = —1, 0, 1 we can avoid dealing 
with f 2 , which however would leave us with the problem of a proof that the solutions are 
integrals. Alternatively, we can choose the range i — 0, 1,2, and then the solutions are 
automatically integrals, as soon as it is established that they are even functions of e. 

A symbolic solution of the system consisting of 18 linear equations, namely eqs. (40)- 
(42) with % — 0, 1,2 along with nine simple equations (56 )— (61 ), would require astronom- 
ical amounts of memory, because of the complexity of f 2 . However, this task becomes 
manageable and even simple for fixed (numerical) values of the phase variables (m,p) and 
of the parameters Wj, while leaving e a symbolic variable. For rational values of mfc,pfc, ujf, 
all computations can be done precisely (in rational arithmetic). This means that atj, ftj, 
and 7j can be evaluated, as functions of e, at arbitrary points in Q 9 (m,p, ui, u 2 , uj 3 ). A 
big number of such evaluations provides us with a convincing evidence in favor of the 
claim that these functions are even in e. 

In order to obtain a rigorous proof without dealing with f 2 , further linear relations 
would be necessary. Before introducing these, we present some preliminary considera- 
tions. Assuming that $i,$2 ; ( ^ ) 3 are HK-bases with one- dimensional null-spaces, results 
of Theorem 13 on the HK-basis ^ tell us that the row vector (0^7, ds, dg) is the unique left 



null- vector for the matrix 



Mo 
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normalized so that 



(d 7 , 4, d 9 )Mi = (1,1,1), where M 1 



OC\ a>2 «3 

fa fa fa 
v7l 72 73, 



Note that due to eqs. (56)-(59) the matrix M 2 has at most four (linearly) independent 
entries. Denoting the common values in these equations by A, B, C, D, respectively, we 
find: 





( "4 


«5 


a 6 \ 




M 2 = | 


fa 


fa 


fa 


H 




V74 


75 


76 J 





D A/fa — Us) A/(ujx — lo 2 ) 
j 2 -uj 3 ) L> B/(uj 2 -ujx) 

Jl-LUo) C/(u)x—wi) D 



(62) 



The existence of the left null-vector (cf 7 , g? 8 , <ig) shows that det(M 2 ) = 0, or, equivalently, 
AB BC CA 



(^1-^3)^2-^3) (ui — ux){uz — u\) (uj 3 - u 2 ) (coi - to 2 ) 



0. 



(63) 



From eqs. (62) and (63) one easily derives that the row 

A 



B 



C 



C 



D , D — 

UJ2 — OJ3 UJ3 — UJ 2 U)\ — Co> 3 UJ3 — UJ\ 

= (a 4 -fa- 74, -as +fa~ 75, -«6 ~ fa + 7e) 



D 



A 



B 



UJ\ — Ld 2 UJ 2 — UJ\ 



is a left null-vector of the matrix M 2 , and therefore (d-r,ds,dg) is proportional to this 
vector. The proportionality coefficient can be now determined with the help of the PSLQ 
algorithm and turns out to be extremely simple. Namely, the following relations hold: 



«4 - fa ~ 74 
-OC5 + fa ~ 75 

-a 6 - fa + 7 6 



D 
D 
D 



B 


-C 


1 


u) 2 


- ^3 


2 


C 


-A 


1 


^3 


- UJx 


2 


A 


- B 


1 


u x 


- UJ 2 


2 



d 7 , 
d 9 . 



(64) 
(65) 
(66) 



Only two of them are independent, because of eq. (46). We note also that, according to 
eq. (53 ), one has 



«4 + fa + 74 
"5 + fa + 75 
"6 + fa + 76 



D + 
D + 
D + 



B-C 



UJ 2 


- ^3 


c 


-A 


^3 


- UJx 


A 


- B 


UJx 


— U) 2 



9 4, 
95, 
96- 



(67) 
(68) 
(69) 



Equations (64)-(69) and (63) are already enough to determine all four integrals A, B, C, D, 
that is, all ctj, fa, 7^ with j = 4, 5, 6, provided it is proven that they are indeed integrals. 
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These (conditional) results read: 

1 + e 2^(2) + e 4 A (4) + e 6 A (6) + e S A (S) 

A -- 
B -- 
C 

D -- 



2e 2 A 






1 + e 2 B {2) + ,4^(4) + e i 


5(6) + 


5(8) 


2e 2 A 






1 + e 2C<(2) + ,4^(4) + £ 6 


C(6) + e 8 


C(8) 


2e 2 A 






p\+pI+pI + e 2 ^ (4) + 


e 4 D (6) + 


e6£)(8) 



2A 



(70) 
(71) 
(72) 
(73) 



where A^ 2q \ B^ 2q \ C^ 2q \ 5^ are homogeneous polynomials of degree 2q in phase vari- 
ables, for instance, 



A (2) 
D (4) 



5(2) = C {2) 

m\ + ml + ml + (uj 2 + ^ 3 - 2uj{)p\ + (u; 3 + u) x — 2uj 2 )p\ + (coi + uj 2 - 2lu 3 )pI, 
(mipi + m 2 p2 + m 3 p 3 ) 2 

+ {p\ +Pl+ Pi) ((^2 + ^3 - 2ui)pl + (w 3 + wi - 2cu 2 )P2 + (wi + uj 2 - 2^3)^3) • 



We remark that eq. (63) tells us that no more than three of the functions A,B,C,D 
are actually functionally independent. Computation with gradients shows that A,B,C 
are functionally independent, indeed. Moreover, all other previously found integrals J, 
d,7,ds,dQ, and gi,g 2 ,g3 are functionally dependent on these ones. 



Theorem 16 The sets (34)-(36) are HK-bases for the map (31) with dim K^, 1 (m,p) 
dim K<s, 2 (m,p) = dim K^ 3 (m,p) = 1. At each point (m,p) G M 6 there holds: 



K® 2 (m,p) 
K$ 3 (m,p) 



[«i 
[Pi 
[7i 



a 2 : a 3 : 0J4 : 0J5 : a§ : — 1], 
72 : 73 : 74 : 75 : 76 : -1], 



where aj,/3j, and jj are rational functions of (m,p), even with respect to e. T/iey are 
integrals of motion for the map (31 ) and satisfy linear relations ( 56)-(61 ). For j = 4, 5, 6 
they are given by eqs. (62), (72), (75). For j = 1,2,3 they are of the form 

h (2) + e 2 h (4) + ^(6) + e 6 h (8) + e8/l (10) + e 10^(12) 



2e 2 (p 2 +p 2 +p 2 )A 



(74) 



where h stands for any of the functions aj,[3j,jj, j = 1,2,3, and the corresponding h^ 2q} 
are homogeneous polynomials in phase variables of degree 2q. For instance, 



(75) 



«! 2) = 




a? 


= -h, 


(2) 

a 3 ' 


= -h, 


A 2) = 


-h, 


& 


= H 3 — I 2 , 




= -h, 


7i W = 


-h, 


7f 


= -h, 




= H 3 -I 3 



The four functions J, ct\, (3\ and 71 are functionally independent. 
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Proof. The proof consists of several steps. 

Step 1. Consider the system for 18 unknowns ctj, /3j,7 3 -, j = 1,...,6, consisting of 



17 linear equations: eqs. (40)-(42) with i — 0, 1, eqs. (56)-(61), and eqs. (64), (65). 
This system is underdetermined, so that in principle it admits a one-parameter family of 
solutions. Remarkably, the symbolic MAPLE solution shows that all variables a 3 -,/3j,7j 
with j = 4,5,6 are determined by this system uniquely, the results coinciding with eqs. 



(62), (72)-(73). (Actually, the MAPLE answers are much more complicated, and their 
simplification has been performed with SINGULAR, which was used to cancel out com- 
mon factors from the huge expressions in numerators and denominators of these rational 
functions.) Since these uniquely determined aj,(3j,jj with j = 4,5,6 are even functions 
of e, this proves that they (i.e., A, B, C, D) are integrals of motion. 

Step 2. Having determined atj,[3j,jj with j = 4,5,6, we are in a position to compute 
atj,Pj,jj with j = 1,2,3. For instance, to obtain the values of cnj with j = 1,2,3, we 



consider the symmetric linear system (40) with i = —1,0,1 (and with already found 
ot±, a*,, a^) . This system has been solved by MAPLE. The solutions are huge rational 
functions which however turn out to admit massive cancellations. These cancellations 
have been performed with the help of SINGULAR. The resulting expressions for a\, a 2 , «3 



turn out to satisfy the ansatz (74) with the leading terms given in the first line of eq. 
(75). (All further terms can be found in [Worksheets].) However, this computation does 
not prove that the functions so obtained are indeed integrals of motion. To prove this, 
one could, in principle, either check directly the identities ctj o f — atj, j — 1, 2, 3, or verify 



equation (40) with i = 2. Both ways are prohibitively expensive, so that we have to look 
for an alternative one. 

Step 3. The results of Step 2 yield an explicit expression for the function 

F — (u 2 — oo 3 )ai + (uj 3 - uji)a 2 + {uj\ - ^2)0:3, (76) 

which is of the form 

(tg - cu 3 )(l + t 2 F^ + e 4 F^ + e 6 F( 6 ) + e 8 F^) 

It is of a crucial importance for our purposes that it can be proven directly that F is an 
integral of motion. We have proved this with the method (G) based on the Grobner basis 
for the ideal generated by discrete equations of motion. The application of this method 
to F is more feasible that to any single of atj, j = 1,2,3, because of the cancellation of 
the huge polynomial coefficient of e 10 in the numerator of F. Actually, more is true: F 
is not only an integral, but is functionally dependent on the previously found ones, say 
on J,d 7 ,d 8 . For a proof of this claim, it would be most favorable to find the explicit 
dependence F = F(J, d?, d$), but it remains unknown to us. Instead, we have chosen the 
way of verification that 

VF G span(VJ, Vd 7 , Vd 8 ). 

This is easily checked numerically for arbitrarily many (rational) values of the data in- 
volved. For a symbolic check, one has to prove the existence of three scalar functions 
Ai, A2, A 3 such that 

VF = Ai V J + A 2 Vd 7 + A 3 Vc/ 8 - 
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This is the system of six equations for three unknowns. Since J does not depend on m&, 
one can determine A2, A3 from a system of only three equations: 



V m -F — \ 2 W m d7 + A3 V m dg. 



After that, 



it remains to check that V P F 



\ 2 V p d 7 



\ 3 V p d 8 is proportional to V P J. 



Clearly, these computations can be arranged so as to verify vanishing of certain (very 
big) polynomials. We have been able to perform these computations with the help of 
SINGULAR for symbolic rrik,Pk but with (several sets of) numeric values of coefficients 
u) k only. 

Step 4- The result of Step 3 allows us to proceed as follows. Consider the system of 
three linear equations for oil, a 2 , a 3 , consisting of (40) with i = 0,1, and of 



(u 2 — uJ-ijdi + (cc>3 — coi)a 2 + {0J1 — uj 2 )a 3 = F, 

where F is the explicit expression obtained and proven to be an integral on Step 3. This 
system can now be solved by MAPLE; the results, again simplified with SINGULAR, 
are even functions of e (actually, the same ones obtained on Step 1 from the symmetric 
system). Non-even polynomials in e of degree 7 cancel in a miraculous way from the 



numerators and the denominator. Now Proposition [14] assures that these solutions are 
integrals of motion. 

Step 5. Finally, in order to find (3\, (5 2 , P3 and 71,72,73, we solve the two systems 
consisting of (41), resp. (42) with i — 0,1, and the first, resp. the second linear relation 
in eq. (60). The results are even functions of e, satisfying the ansatz (74) with the leading 



terms given in eq. (75). Proposition 14 yields that also these functions are integrals of 
motion. □ 
MAPLE worksheets for all computations used in this section can be found in [Work- 
sheets] . 



6.9 Preliminary results on the Hirota-Kimura-type discretiza- 
tion of the general flow of the Clebsch system 

The general flow of the Clebsch system, depending on three real parameters b\, b 2 , b 3 (or, 
rather, on their differences fej — bj, which gives two independent real parameters), reads 
as follows: 

m = m x Cm + p x Bp , 

( 77 ) 

p = p x Cm, 

where B = diag(&x, b 2 , b 3 ) and C = diag(ci, c 2 , c 3 ) with 

b 2 -63 63 - 61 h - b 2 
ci = , c 2 = , c 3 = . (78) 

LU 2 — CU3 LU3 — UJ\ LUi — 0J 2 

This flow is Hamiltonian with the quadratic Hamilton function 

H= l -(m, Cm) + l - (p, Bp) = l - ^{c k m 2 k + b k p 2 k ) = ^(hh + b 2 I 2 + 63/3). 

fc=i 
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In components eqs. (77) read: 



rrii 
m 2 
m 3 
Pi 

h 
r>3 



(c 3 - c 2 )m 2 m 3 + (6 3 - b 2 )p 2 p 3 , 
(ci - c 3 )m 3 m 1 + (61 - b 3 )p 3 pi, 
(c 2 - ci)m l m 2 + (6 2 - &i)piP2- 
c 3 m 3 p 2 - c 2 m 2 p 3 , 
c x m\p 3 - c 3 m 3 pi, 
c 2 m 2 pi - c x mxp 2 . 



The KH discretization of the flow (77) reads 



m — m 



p — p 



In components: 



mi — mi 
wi2 — m 2 
m 3 — m 3 

pi ~Pi 

P2 ~Pl 
P3-P3 



e(m x Cm + m x Cm + px Bp + p x Bp 
e (p x Cm +px Cm) . 



e(c 3 - c 2 )(m 2 m 3 + m 2 m 3 ) + e(6 3 - b 2 ) (p 2 p 3 + p 2 p 3 ) , 
e(ci - c 3 )(m 3 mi + m 3 m x ) + e(b x - b 3 )(p 3 pi +p 3 pi), 
e(c 2 - ci){m l m 2 + mim 2 ) + e(b 2 - &i)(pip 2 +P1P2), 
ec 3 {fh 3 p 2 + m 3 p 2 ) - ec 2 (m 2 p 3 + m 2 p 3 ), 
eci(mip 3 + mip 3 ) - ec 3 (m 3 pi + m 3 pi), 
ec 2 (m 2 pi + m 2 p~i) - eci{mip 2 + m^)- 



(79) 



In what follows, we will use the abbreviations i>j 



bj and Cjj = q — Cj. The linear 



system (79) defines an explicit, birational map / : 



f(m,p,e) 



M~ 



where 



M(m,p, e) 



/ 1 ec 23 m 3 ec 23 m 2 

ec 31 m 3 1 ec 3 im! 

ec i2 m 2 ec i2 mi 1 

ec 2 p 3 -ec 3 p 2 

-ecip 3 ec 3 pi 

\ ecip 2 -ec 2 pi 



(m,p,e) 





eb 31 p 3 
tb 12 p 2 
1 

ec 3 m 3 
-ec 2 m 2 



tb 23 p 3 


eb 12 Pi 
-ec 3 m 3 
1 

ecimi 



10) 



eb 23 p 2 \ 
eb 31 pi 


ec 2 m 2 
— ecimi 

1 / 



As usual, map (80) possesses the reversibility property 

/ _1 (w,P,e) = f(m,p, -e). 



Conjecture 17 All claims of Theorems 10, 11 hold also for the discretization (80) of the 

(81) 



general flow of the Clebsch system, with the HK-basis $0 being replaced by 



$0 = (Pi,P 2 ,P 3 ,TOi,m 2 ,m 3 ,l) 



33 



This conjecture is supported by numerical results based on the algorithm (N). The claim 



concerning the set <£>o given in eq. (81) is proven symbolically. In order to keep the 



notations compact, we give here this proof for the second flow of the Clebsch system only. 

Recall that the first flow of the Clebsch system, considered in Sect. [6j corresponds to 
hi = LOi and c$ = 1. The second flow is characterized by the Hamilton function 

H = ]-H 2 = ^(uJ%ml + u 2 m 2 2 + u 3 m 2 3 - u 2 u 3 pj - uJiU 3 p 2 2 - UiU 2 p 2 3 ). 
In other words, the choice of parameters bk characterizing the second flow is 

h = -u} 2 u 3 , b 2 = -o> 3 o>i, 6 3 = -Wio> 2 , (82) 

so that 

ci = u>i, c 2 = uj 2 , c 3 = uj 3 . (83) 

For the HK discretization of the second Clebsch flow, we give a more concrete formulation 
of our findings concerning the HK-basis $o> including a "nice" integral. 



Theorem 18 For the map (80) the set of functions (81) is a HK-basis with dim K^ (m,p) - 
1. At each point (m,p) e M. 6 there holds: 

K<s> (m,p) = [ei : e 2 : e 3 : e 4 : e 5 : e 6 : -1], 
where all are fractional-linear functions of a single integral L = L(m,p,e) of the map 



(80) which is a quotient of two quadratic polynomials in rrik,Pk- 



If the coefficients &fc,c& are as in eqs. (8ty ), (83), then the integral L can be taken as 
Ei{uJim\ + u 2 uj 3 p\) + E 2 [uj 2 m\ + u 3 Uijp$) + E 3 (oj 3 m 3 + uj^pf) 



1 + e 2 uji(ujim\ + u 2 uj 3 pf) + e 2 uj 2 (u 2 ml + u 3 Uipf) + e 2 uj 3 {uj 3 ml + ou^pf) 



2\ ' 



with 



Ei = uj 3 uJi + UJ\UJ 2 — uj 2 lu 3 , E 2 = ujiuj 2 + uj 2 uj 3 — UJ3UJ1, E 3 = uj 2 uj 3 + uj 3 lu 1 — cu 1 lu 2 . (84) 

In this case there holds: 
ei e 4 



uj 2 uj 3 uj\ 



e 2 e 5 



Ei 2 

Ei + e 2 (uJi - u 2 )E 3 (uj 2 m 2 2 + u 3 Uip 2 2 ) + e 2 {ui - u 3 )E 2 {u 3 ml + Uiuj 2 p\) 



Ei(ujim\ + 0J 2 uj 3 p 2 i) + E 2 (uj 2 ml + uj 3 ujipI) + E 3 (uj 3 ml + 0JiUJ 2 p 2 3 ) 



2\ ' 



E- 



2 .2 



UJ 3 UJl L0 2 



e 3 e 6 



e lo 2 



E 2 + e 2 (u) 2 - uj 3 )Ei(uj 3 ml + uoiuo 2 pl) + e 2 (tu 2 - u)i)E 3 (uiim\ + uo 2 uo 3 p\] 
Ei(uim\ + uj 2 uj 3 p\) + ^2(^2^2 + ^s^ipf) + E 3 (uo 3 ml + w^pl) 



2\ ' 



E- 



UJlU) 2 UJ 3 



3 2 

e u; 3 



E 3 + e 2 (io 3 - uJi)E 2 (ujim 2 i + 0J 2 uJ 3 p\) + e 2 [uj 3 - uj 2 )Ei{uj 2 m\ + uJ 3 uJip 2 2 ) 
E x {uim{ + u^Pi) + E 2 {y) 2 m\ + o; 3 ^iPi) + £3(^3™! + uiu 2 p 2 3 ) 
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The numerator of L, 



L(m,p, 0) = Ei(uJim\ + uj 2 uj 3 p\) + E 2 (uj 2 m 2 2 + uj 3 ujip\) + E 3 (uj 3 m\ + ujiuj 2 pl), 

is a linear combination of quadratic integrals of motion of the continuous Clebsch system. 

Proof. We will only present the proof for the second Clebsch flow. The claim of the 
theorem refers to the linear system 



{exp\ + e 2 p\ + e 3 p 2 3 + e±m\ + e 5 m 2 2 + e 6 m 2 3 ) o f l (m,p) = 



1 



for i from the ranges containing 6 consecutive numbers, such as i G [—2, 3] or i G [—3, 2]. 
As the solution of such a system clearly requires more iterates of the map / than could 
be handled symbolically, we follow recipe (E) and look for linear relations between ej. It 
turns out to be possible to identify the following five relations: 

uiiei — uj 2 ijj 3 e^ = 0, (85) 

uj 2 e 2 - oj 3 uJie 5 = 0, (86) 

u 3 e 3 - UxUJ 2 e & = 0, (87) 

d - e 2 - (ui - uo 2 )e G = e 2 Lul(uJi - oj 2 ), (88) 

e 3 - e x - (u 3 - uJi)e 5 = e 2 ul{u 3 - wi). (89) 

Of course, there holds also a third non-homogeneous relation: 

e 2 - e 3 - (u 2 - uj 3 )e 4 = e 2 u 2 (u 2 - io 3 ), 

but actually it is a consequence of the previous five. As usual, these (at this point 
conjectural) identities can be (and have been) found using the PSLQ algorithm. Now 
we obtain the six functions by solving a simple system of six linear equations which 
involves no iterates of the map / at all and consists of 

eiPx + e 2P 2 + CsPi + 64^1 + e 5 ml + e§m\ = 1 



along with the relations (85 )-(89 ). The solution is given in the formulation of the theorem. 
To prove that the function L is an integral of motion one can use a straightforward 
computation using MAPLE. Also a proof based on the equations of motion alone can 



be given, similar to the proof for L (see proof of Theorem 12). The last claim of the 
theorem about L(x, 0) follows in the limit e — > 0, but can be also easily checked directly, 
by verifying conditions (78) for 6« = UjUkEj and q = uJiEi with E{ from eq. (84). These 



conditions are satisfied due to the identities 

ujjEi — cjiEj = (u>i — LOj)Ef., 
where k) is any permutation of (1,2,3). □ 
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7 Conclusions 



We established the integrability of the Hirota-Kimura-type discretization of the Clebsch 
system, in the sense of 

• existence, for every initial point (m,p) G 1R 6 , of a four-dimensional pencil of quadrics 
containing the orbit of this point; in our terminology, this can be formulated as 
existence of a HK-basis with a four-dimensional null-space, consisting of quadratic 
monomials; 

• existence of four functionally independent integrals of motion (conserved quantities). 

Numerical experiments show that this remains true also for an arbitrary flow of the 
Clebsch system. It is interesting to remark that the maps generated by Hirota-Kimura 
discretizations of various flows do not commute with each other. It would be important 
to understand whether some analog of commutativity of the continuous flows survives in 
the discrete situation. 

Our investigations were based mainly on computer experiments. Our proofs are com- 
puter assisted and were obtained with the help of symbolic calculations with MAPLE, 
SINGULAR and FORM. A general structure behind these facts, which would provide 
us with more systematic and less computational proofs and with more insight, remains 
unknown. In particular, nothing like a Lax representation has been found. Nothing is 
known about the existence of an invariant Poisson structure for these maps. (For a sim- 
pler system, Hirota-Kimura discretization of the Euler top, an invariant volume measure 
as well as a bi-Hamiltonian structure have been found in [Petrera and Suris 2008].) 

Hirota and Kimura demonstrated that their discretization leads to an integrable map 
also for the Lagrange top [Kimura and Hirota 2000]. Our preliminary investigations 
have shown remarkable features pointing towards the integrability of the Hirota-Kimura 
discretizations of the following systems: Zhukovsky-Volterra gyrostat; so(4) Euler top 
and its commuting flows; Volterra and Toda lattices; classical Gaudin magnet. Based on 
these observations, we formulate the following hypothesis. 

Conjecture 19 For any algebraically completely integrable system with a quadratic vector 
field, its Hirota-Kimura discretization remains algebraically completely integrable. 

If true, this statement could be related to addition theorems for multi-dimensional 
theta-functions. Such a relation has been already established for the Hirota-Kimura dis- 
cretization of the Euler top, which can be solved explicitly in elliptic functions [Suris 2008]. 
In our ongoing investigations, we hope to establish integrability of the above mentioned 
discrete time systems and to uncover general mechanisms behind it. 
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